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Abstract 

Even though our theoretical understanding of dendritic solidification is rela- 
tively well developed, our current ability to model this process quantitatively 
remains extremely limited. This is due to the fact that the morphological 
development of dendrites depends sensitively on the degree of anisotropy of 
capillary and/or kinetic properties of the solid- liquid interface, which is not 
precisely known for materials of metallurgical interest. Here we simulate the 
crystallization of highly undercooled nickel melts using a computationally effi- 
cient phase-field model together with anisotropic properties recently predicted 
by molecular dynamics simulations. The results are compared to experimen- 
tal data and to the predictions of a linearized solvability theory that includes 
both capillary and kinetic effects at the interface. 

*Present address: Department of Radiation Oncology, Yonsei University, College of Medicine, Seoul 
132, Korea. 
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I. INTRODUCTION AND OVERVIEW 



The freezing of a supercooled liquid or a supersaturated liquid mixture of two or more 
components generally leads to the formation of highly branched dendritic crystals in materi- 
als where the solid-liquid interface is rough on an atomic scale, but smooth on a macroscopic 
scale M . This process is of considerable practical importance to metallurgists because com- 
mercial metallic alloys often form dendrites during freezing, and the resulting microstructure 
controls the properties of the final solidification product . Moreover, dendritic growth has 
been of fundamental interest to physicists in the general context of pattern formation in 
non-equilibrium dissipative systems || [| . 

The solidification of a pure undercooled melt is one of the simplest methods to form den- 
dritic crystals and to investigate their growth behavior and morphology under well-controlled 
conditions. For this reason, this process has been extensively studied both theoretically and 
experimentally for several decades. Traditionally, experimental studies have focused on the 
characterization of the unique relationship between the steady-state growth velocity V of 
the dendrite tip and the bulk undercooling AT = Tm — T^, where Tm is the melting point 
temperature and is the initial temperature of the melt. 

A first major step in the development of the theory of dendritic growth was made over 
half a century ago by Ivantsov || who analyzed this problem in the simplifying limit where 
both capillary and kinetic effects at the interface are neglected. In this approximation, the 
solid-liquid interface grows at the melting temperature and the only rate-limiting process is 
the diffusive transport of the latent heat of freezing into the undercooled liquid. The basic 
equations governing the growth are then the diffusion equation 

and the condition of heat conservation at the growing interface 

LV n = c p Dh-(yT\ - VT| ), (2) 
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where V n is the normal velocity of the interface, h is the directional normal to the interface 
pointing into the liquid and n ■ VT\s,l is the normal gradient of temperature on the solid 
(S) or liquid (L) side of the interface; also, L and c p are, respectively, the latent heat of 
freezing and the specific heat at constant pressure per unit volume, and D is the thermal 
diffusivity (with both c p and D taken here to be the same in both phases). 

Ivantsov looked for needle-shape steady-state solutions of Eqs. [I] and || growing at a 
constant velocity V. The assumption that growth takes place in steady-state allows one 
to replace dT/dt in Eq. [l] by -VdT/d 2 in a frame of reference moving with respect to 
the sample at speed V, where +z is the growth direction. The constraint that the needle 
shape is preserved in time implies that V n = V^cos^, where 9 is the angle between n and 
the z-axis. Remarkably, Ivantsov found that paraboloids of revolution are exact solutions 
to this steady-state growth problem for an isothermal interface at T = Tm- Moreover, he 
derived a relationship between the Peclet number, P = pV/(2D), where p is the dendrite 
tip radius, and the dimensionless undercooling AT/(L/c p ). Three decades later, Glicksman 
et al. [7| demonstrated by accurate measurements in a transparent organic system that this 
relation is well satisfied. 

The Ivantsov relation does not predict V and p independently, but only their product. 
Consequently, theoretical studies subsequent to Ivantsov's work concentrated on understand- 
ing how this velocity is uniquely determined by including capillary and kinetic effects. These 
effects make the interface non-isothermal, with a temperature determined by the classic 
Gibbs-Thomson condition augmented to allow for a departure from local equilibrium 



Tj = T M -^fJ2 

L i=l,2 
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The second term on the right-hand-side of Eq. [3] represents the purely equilibrium effect 
of capillary forces on the temperature of a curved interface, where 7(n) is the excess free- 
energy of the solid-liquid interface, 9i are the local angles between the normal direction 
n and the two local principal directions on the interface, and Ri are the principal radii 
of curvature. This term has the crucial effect that it makes the interface stable against 



short-wavelength perturbations. The last term on the right-hand-side of Eq. |] is the non- 
equilibrium undercooling of the interface that drives crystallization by attachment of atoms 
at the interface, where fi(n) is a kinetic coefficient. The solution of the steady-state growth 
problem becomes highly non-trivial when the isothermal condition T = Tm is replaced by 
Eq. ^, and several dendrite growth theories after Ivantsov (reviewed in Ref. made some 
uncontrolled approximations that were later proven to be incorrect. 

The fundamentally correct solution to this problem emerged in the early nineteen eight- 
ies from the study of simplified models where the interface dynamics is governed by local 
growth laws (a purely geometrical model where V n depends on the curvature and its second 
derivative with respect to arclength along the interface ||, and a boundary-layer model || 
that is an approximation of the dendritic growth equations in the limit where the thickness 
of the thermal diffusion layer surrounding the tip is small compared to p). These models 
had the advantage that they could be simulated numerically and studied analytically with- 
out uncontrolled approximations that had been made in previous theoretical attempts to 
tackle directly the dendrite growth equations with the full non-locality of the diffusion field. 
Building on insights gained from these studies, solutions to the the full equations were then 
subsequently obtained in two [Tt]-[rj| and three dimensions PT5| , [TB . 



The key insight gained from these studies is that dendritic evolution is controlled by 
a delicate balance between microscopic capillary effects and macroscopic heat transport in 
which the anisotropic properties of the solid-liquid interface play a crucial, and mathemati- 
cally subtle, role. First, capillary forces act as a singular perturbation that transforms the 
Ivantsov continuous family of needle crystal solutions into a discrete set of solutions. Second, 
capillary or/and kinetic effects must be anisotropic for solutions to exist; in this case only 



the solution with the largest velocity in this discrete set is linearly stable [|1y[] and therefore 
a potential candidate for the observed dendrite tip. 

The self-consistent calculation used to determine this velocity involves a solvability con- 
dition, which is simply the condition that a physically admissible solution to the steady-state 
growth equations must be smooth at the tip. Consequently, the term solvability theory is 



commonly used to refer to this type of calculation whether it is carried out numerically using 
a boundary integral approach [|IT|- |i3| , |H)[ | or analytically |IO| , [i^ , |i6|| using a WKB (Wentzel- 



Kramers-Brillouin) approximation in the limit of weak anisotropy. 

Over the last decade, new powerful computational approaches have been developed to 
simulate the full time-dependent evolution of the solid-liquid interface governed by Eqs. [I|-|3|. 
In particular, the phase-field method [|l8j has emerged as a powerful algorithm to simulate 



this evolution in both two and three dimensions ||19|-Pq] . This approach has the well-known 
advantage that it facilitates front tracking by making the phase boundaries spatially diffuse 
with the help of scalar fields that naturally distinguish between different phases. Using 
this approach, it has been possible to simulate the morphological development of dendrites 
and to validate quantitatively that the dynamically selected dendrite tip operating state is 
indeed the one predicted by solvability theory |22] . 



Despite this progress, our current ability to model dendritic growth quantitatively for a 
given material remains very limited. This is mainly due to the lack of knowledge of capillary 
and kinetic properties (7(1^) and /x(n), respectively, which enter through Eq. |3|) and, in 
particular, of their anisotropies that strongly influence the interface dynamics. Exceptions 
are some transparent organic materials P?^pP| and only one metal system (Al-Cu) |SD| for 
which values of the capillary anisotropy have been estimated from experimental measure- 
ments of equilibrium shapes. Even for these materials, however, interface kinetic effects have 
not been quantified even though they could be important. Traditionally, kinetic effects have 
been assumed to be negligibly small compared to capillary effects at small undercoolings, 
where the growth rate is small, and the opposite for large undercoolings. However, to predict 
dendritic evolution over the entire experimentally accessible undercooling range requires a 
precise knowledge of both capillary and kinetic properties. 

The hope to model dendritic evolution quantitatively in metal systems has been recently 
revived by progress in modeling interfacial properties using large scale molecular dynamics 
(MD) simulations pH~53|. The interface kinetic coefficient has been calculated for Ni and Cu 



by crystallizing large slabs of liquid with interatomic potentials derived from the embedded 



atom method |3T| . Values of this coefficient were found to be in good quantitative agreement 



with a model of Broughton et al. |35|] used previously to interpret simulations with the 
Lennard- Jones potential, and about half an order of magnitude smaller than the value 
predicted by a model of Turnbull [[36] that has been used in the metallurgical literature 
to model rapid dendritic solidification. Furthermore, the interface kinetics was found to be 
strongly anisotropic, with growth at the same interface undercooling being faster along [100] 
than [110], and faster along [110] than [111]. A similar magnitude of anisotropy has also 
been found in simulations with the Lennard- Jones potential [[37] . 



The anisotropy of the interfacial energy is more difficult to compute for materials with 
rough interfaces because it is extremely weak, the maximum variation of 7(n) over all orien- 
tations (n) of the interface being typically of the order of one percent. Very recently, however, 
it has been demonstrated that this anisotropy can be accurately computed by monitoring 
interfacial fluctuations during MD simulations and extracting the interfacial stiffness, i.e. 7 
plus its second derivative with respect to orientation in a given plane, which is one order of 
magnitude more anisotropic than 7 itself 

The goal of the present work is to model dendritic growth in undercooled metallic melts 
fully quantitatively by linking atomistic and phase-field simulations. We modify the existing 
phase-field formulation of the crystallization of a pure melt in order to achieve an efficient and 
accurate modeling of the large-undercooling regime where interface kinetic effects become 



important, and we use a recently developed multi-scale algorithm [[26] to simulate this model 
in three dimensions. Furthermore, we use as input into this model the anisotropic capillary 
and kinetic properties of the solid-liquid interface that have been recently predicted by 
MD simulations [j31|jn^ . This allows us to make quantitative predictions of morphological 
development that can be directly compared with experiments. 

We focus on pure Ni for which the interface properties have been most accurately modeled 
to date by MD simulation. This system also has the advantage of having been extensively 
studied experimentally. Various groups have measured the relationship between the growth 
velocity of the interface and the bulk undercooling [j38j-fl2l . Moreover, the morphological evo- 
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lution of the envelope of the solidification front during recalescence (rapid freezing following 
triggered nucleation) has been investigated by thermal imaging pT| , [42|| . 

We write down the phase-field model and examine its sharp interface limit in section II. 
We then discuss the numerical implementation of this model in section III. Next, we discuss 
the results of phase-field simulations in section IV. We mainly focus here on the comparison 
of these results with experimental observations. The comparison with the predictions of 
solvability theory is only briefly summarized and will be presented in more detail elsewhere. 
Finally, conclusions and prospects are given in section V. 



II. PHASE-FIELD MODEL 



A. Basic equations 



Over the last decade, the phase-field method [ffH] has emerged as a powerful computa- 



tional approach to simulate morphological development during solidification [|i9| -|26ll. This 
method has the well-known advantage that it avoids to explicitly track a sharp phase bound- 
ary by making the interface region between solid and liquid spatially diffuse over some finite 
thickness. The general form of the phase-field model that has been widely used to simulate 
the solidification of a pure melt is defined by the equations 

, 86 5F 



r(n) 



dt 
du 
~dt 



(4) 
(5) 



where 8F/ 5(f) denotes the functional derivative with respect to (f> of the free energy of the 
two-phase system. The latter assumes the phenomenological form 

~W 2 (n) 



dx 



■|V0| 2 + /^(0, M ) 



(6) 



where the integral in Eq. ||] is over the volume of the system and we have defined the 
dimensionless temperature field 
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u^{T-T M )/{L/c p ). (7) 

Here, fd w is dimensionless and W has dimension of length; consequently, F has dimension 
of volume. It is important to emphasize that there is a significant amount of freedom in the 
choice of the functional form of the free energy density fd w {<f),u) in the phase-field model. 
The first main requirement is that fdw should have the form of a double well potential with 
two minima corresponding to the liquid and solid phases, which are chosen here at <fi = ±1 
with the plus (minus) sign corresponding to the solid (liquid). Thermodynamically, this 
guarantees the existence of a well-defined solid-liquid interface of thickness ~ W, i.e. where 
varies between +1 and —1 over this thickness along a direction that is locally normal to 
the interface. The second physical requirement is that the difference in bulk free-energy 
between solid and liquid, which drives the phase transformation, should be a monotonously 
increasing function of the interface undercooling; this is equivalent, here, to requiring that 
fdw(+^,u) — fdw(—^,u) be a monotonously increasing function of u. A third requirement is 
that the two minima of fdw remain at ±1 for different temperatures above and below the 
melting point. A choice that satisfies all three requirements is |pll- p6i^3 



fdw(<P,u) = f (</>) + \ug(<f>), (8) 

with A a dimensionless constant, /(</>) = — 4> 2 /2 + 4 /4, and 

g(<f>) =4>- 20 3 /3 + 5 /5. (9) 

This form has been widely used to model dendritic solidification quantitatively in a low 
undercooling regime where interface kinetic effects are vanishingly small [^,^,^]. Here, 
however, we focus on the opposite high undercooling regime where kinetic effects are domi- 
nant. In order to model accurately the kinetics of the interface in this regime, it turns out 
to be necessary to use the more general form 

f dw (<f>,u) = f (<!>)+ h(\u)g(<i>), (10) 

which still satisfies the above requirements if h is a monotonously increasing function of its 
argument. Even on today's computers and with parallel codes, phase-field computations 
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are only feasible for values of W much larger than the real atomic-scale width of the solid- 
liquid interface, denoted here by 5. With the double well potential defined by Eq. [8] and 
W 3> 5, the velocity of a planar interface in the phase-field model is a nonlinearly increasing 
function of the interface undercooling for the high velocity range (30 — 60 m/sec) that has 
been measured experimentally in Ni |33-H2[. In contrast, atomistic simulations of the same 



material |yj have shown that the relationship between the growth rate and the interface 
undercooling is still linear in this velocity range. This problem can be solved by introducing 
a function h in the double well potential defined by Eq. |10| that allows one to freely choose 
the relation between the interface undercooling and the thermodynamic driving force for 
the phase transformation. Hence, even when the relationship between the interface velocity 
and the driving force is nonlinear, one can model a linear dependence of the velocity on 
the interface undercooling. We discuss next how to explicitly compute h and how to choose 
the parameters in the phase-field model by examining the mapping of this model onto the 
sharp-interface equations of the previous section. The incorporation of thermal fluctuations 
into the model is discussed in section III together with the numerical implementation. 



B. Relationship between phase-field and sharp-interface models 

To examine the relationship between the phase-field and the sharp-interface equations 
(Eqs. [l]-|D, it is useful to first rewrite the latter equations in terms of u, which yields 



at = DVu - 



V n = Dn- (Vu| - Vu 



i=l,2 



aJn) + 



d 2 a c (n) 

def 



tL 



(11) 

(12) 
(13) 



where we have defined the microscopic capillary length do = 'JoTmCp/L 2 , which is typically 
on the order of one nanometer, and the kinetic coefficient 



/3(h) = 



(14) 
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Here, 70 is the mean value of the interfacial energy that can be defined by writting j(n) = 
7oa c (n), where the small variation of the interfacial energy with orientation is contained in 
the dimensionless function a c (h) defined in section III below. The interface energy of the 
phase-field model (which has dimension of length because fd w was chosen dimensionless), 

7p /(n) = W(n) f 1 ^2[j{<j>) - /(±1)] d<j>, (15) 



depends only on the shape of the free energy function and is proportional to W . Therefore, 
anisotropic capillary effects are correctly reproduced in the phase-field model by letting the 
interface thickness have the same anisotropy as the interfacial energy ||19|| , that is W(h) = 
Woa c (n) where the direction normal to the interface is defined by n = — V0/|V</»|. 

The phase-field model is well-known to map onto the above set of sharp-interface equa- 
tions in the limit where the interface thickness is small compared to the scale of the mi- 
crostructure. The asymptotic analysis necessary to construct this mapping is rather in- 



volved (see Ref. |22| and earlier references therein). We only mention here the main re- 
sults that are necessary for our computations and to determine h. Let us start with the 
equilibrium Gibbs-Thomson condition. A curved interface in the phase-field model is sta- 
tionary only when the temperature u is constant and the driving force induced by the 
curvature is exactly compensated by the free energy difference between the two bulk phases, 
f dw (+l,u) - f dw (-l,u) = Xu[g(+1) - g{-l)), which yields 



A«kK+i) - g(-i)] = - £ 



ivA n > + 



86? 



;T (16) 



i=l,2 

Using Eq. [15| and comparing the result to Eq. [L3|, we deduce that 

do = ai — , (17) 

where a\ = j\ J2[f(<f>) — /(±l)]/[gr(+l) — g{— 1)] is a constant of order unity (a\ = 5a/2/8 
for the forms of f(4>) and g(<f>) used here). 

Eq. implies that for a given capillary length do, which is the only fixed length scale 
in the dendrite growth problem, one has the freedom to vary the interface thickness Wq 
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in the phase-field model since A is a free parameter. This freedom is due to the fact that 
the free-energy difference on the left-hand-side of Eq. [16] is proportional to A whereas the 
surface energy on the right- hand- side is proportional to Wo- The goal of choosing Wo as 
large as possible in order to make the phase-field equations less stiff is equivalent to choosing 
A ~ Wo /do as large as possible. How large A can be chosen, under the constraint that the 
results be converged (i.e. independent of A), has been shown to depend on the choice of 



assumption made in mapping the phase-field model onto the sharp-interface equations |22 
In the classical asymptotic analysis of the sharp-interface limit of the phase-field model, A is 
assumed to be vanishingly small. This assumption is equivalent, physically, to assuming that 
the temperature is constant in the diffuse interface region in the limit where its thickness 
vanishes. The expression for the kinetic coefficient that results from this analysis is 

A distinct thin-interface limit has also been analyzed where A is assumed to be of order 
unity |]22|| . This analysis computes a correction to the interface kinetics that results from 
the variation of temperature across the interface width, which is assumed finite (or, in simple 
terms, the interface is assumed to be "thin" as opposed to sharp). This analysis yields the 
modified expression for the kinetic coefficient |22| 



. W 2 (n) 
1 - a 2 X- 



Dr in) 



(19) 



where a 2 is a constant of order unity that depends again on the functional forms of the 
functions /(</>) and g(<f>) (a 2 = 0.6267 for the forms used here). This expression has the 
advantage that it remains valid for larger A than the one derived in the classical sharp- 
interface limit. Furthermore, it can even be used to make the interface kinetics vanish by 
choosing r(n) = \a 2 W 2 (h)/D, which makes the square bracket in Eq. |T9| vanish. 

In the early stage of the present investigation, we used the standard form of the double 
well potential given by Eq. R] to simulate dendritic evolution in two dimensions for the 



parameters of pure Ni discussed below, together with Eqs. [17| and [19] to interpret the 
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results. We investigated values of AT/(L/c p ) varying between 0.2 and 0.6, which is the 
undercooling range of interest to make contact with experiments. We found that the dendrite 
tip velocity only became independent of the interface thickness for A of order unity, which 
is equivalent to choosing this thickness comparable to the atomic-scale width 5 of the real 
solid-liquid interface. These simulations turned out to be extremely costly and not feasible 
in three dimensions. The results of the simulations were analyzed in order to elucidate the 
cause of this poor convergence. Plots of the temperature profile along the central axis of the 
dendrite revealed that the magnitude of the temperature variation across the diffuse interface 
is small compared to the magnitude of the total interface undercooling. The latter is itself 
a significant fraction of the total bulk undercooling at this large growth rate. Consequently, 
the correction to the kinetic coefficient predicted by the thin-interface limit is essentially 
negligible. This is easily seen by rewriting Eq. |1^ in the form 

d 



1 — a 2 X- 



(20) 



using Eq. [T7|, where (3q denotes the kinetic coefficient predicted by the sharp-interface limit 
(Eq. |IBD , and where we have neglected anisotropy for simplicity. For the parameters of pure 
Ni predicted by the MD simulations PT| , do/(D(3 ) ~ 1/90, such that the second term inside 



the square bracket in Eq. ^ is small compared to one for A of order unity. 

Further analysis of the results revealed that the linear relationship between interface 
undercooling and velocity, uj = —j3V n , breaks down for the velocity range of interest here. 
Since only the product \u enters the phase-field equation, this breakdown occurs at smaller 
velocity as A is increased. This explains the dependence on A of our dendrite simulation 
results with fd w defined by Eq. ||. To overcome this difficulty, we use the new form of fdw 
defined by Eq. [K]. The velocity of a planar interface is uniquely related to the driving 
force, which is — Xuj in Eq. ^, but —h(\ui) in Eq. [L0| Hence, once we have determined the 
nonlinear relationship between the interface velocity and h, we can then choose h(Xuj) so as 
to recover a linear relationship between interface velocity and uj. Let us first compute the 
nonlinear relationship between the interface velocity and h. This can be done by rewriting 
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Eq. |, with fdw defined by Eq. [10], in a frame translating along the +x direction at velocity 
V n and by making the change of variable y = x/W(n), which yields 

where we have defined the dimensionless interface velocity v = V n r(n)/W(h) and the primes 
on / and g denote differentiation with respect to <p. Physically admissible solutions of Eq. 
|2ll that exhibit a smooth and monotonous variation of <p across the diffuse interface only 
exist for a unique v for a given value of h. They can be found by first writing Eq. [H] as 
a set of two coupled first-order differential equations, and then using a shooting method to 
find the value of v that satisfies the boundary conditions 0(±oo) = =pl. By repeating this 
procedure for different values of h, one can compute the function v(h), which we plot in Fig. 
IJ To make the interface kinetics linear for large A, it now suffices to choose the function 
h(Xui) = v Xui/ai), where v~ l is the inverse of the function v. This inverse function 
is uniquely defined because v is a monotonous function of its argument. Consequently, 
v(h(Xui)) = v (v Xuj/ai)) = —Xui/ai, such that finally 

V n = ^ v(h(X Ul )) = - Ul /f3(h), (22) 
Tin) 



where (3(n) is given by Eq. [18], as desired. For modeling dendrites, the magnitude of Xuj only 
becomes large when the interface is undercooled (uj < 0), whereas this magnitude is small 
during melting (uj > 0), e.g. the remelting of secondary branches. Therefore, we choose 
h(Xui) = v^^—Xui/ai) for uj < and h(Xui) = Xui for m > 0. Finally, by extending the 
above analysis to a moving curved interface, it can also be shown that the Gibbs-Thomson 
condition with a linear kinetic correction term (Eq. [T3| ) is correctly reproduced with the 
new form of fdw defined by Eq. ITOL 



III. NUMERICAL IMPLEMENTATION 



In summary, we simulate the phase-field model defined by Eqs. [| and [5| together with 
the double-well potential defined by Eq. 0, /(<£) = -<p 2 /2 + 4 /4, g{(f) defined by Eq. |, 
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the function h(Xu) defined by h(Xu) = Xu for Xu > and constructed, for Xu < 0, from the 
inverse of the function v(Xu) computed by solving Eq. |21] as outlined above; values of h for 
An < are tabulated beforehand at equally spaced intervals and h is computed by simple 
linear interpolation from these tabulated values during the simulations. The interfacial 
energy and the kinetic coefficient are assumed to have the forms 

4e r 



7(")/7o = a c(n) = (1 - 3e c ) 
P(n)/Po = a k(n) = (1 + 3e fc ) 



1 + 3e fc y 



(23) 
(24) 



consistent with the symmetry of a crystalline material with an underlying cubic symmetry, 
where e c (e^) is the magnitude of the capillary (kinetic) anisotropy. With the above defi- 
nitions, a positive value of e c (e^) favors growth along the [100] directions, which requires 
the interfacial energy (kinetic coefficient) to be maximal (minimal) along these directions. 
Using the definition W(n) = W a c (h) together with Eqs. we obtain 

W(n) = ^a c (n) (25) 

, M 32A 2 d (3o , M /-s / oc x 
T[n) = — a c (n) a k {n) (26) 

that are integrated in the phase- field model using h = — V0/|V0|, or 

r} { r} { r} _ {d<p/dxY + {d<p/d y y + {d<p/d z y 

n x n y n z ' |^| 4 

in Eqs. [23}f?|. 

Thermal fluctuations are incorporated into the model by using the standard Langevin 
formalism that has been widely used in studies of second order phase transitions [44|, and 



more recently in phase- field models of solidification [[25|,[|I|[|6|] . This formalism consists in 
adding random variables to the right-hand side of the phase-field equations, with the noise 
in Eq. |5| written in such a way that the total heat is conserved. These random variables 
are chosen to be uncorrelated in space and time and to obey Gaussian distributions whose 
variances are dictated by the fluctuation-dissipation theorem. This formalism has been used 
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previously by Karma and Rappel |25] to investigate the formation of sidebranches in den- 
dritic solidification, and we follow here exactly the same procedure to include fluctuations. 
The variances of the non-conserved and conserved noises added to Eqs. ^ and [5], respectively, 
are determined by a single dimensionless experimental parameter 

p _ k B Tl[Cp 



where ks is the Boltzmann constant, as detailed in Ref. 25 



Our simulations exploit a hybrid finite-difference diffusion Monte Carlo algorithm that 



integrates efficiently the phase-field equations [pq| . In a first domain D\, consisting of the 
solid plus a thick liquid layer surrounding the solid-liquid interface, the stochastic phase- 
field equations (with the Langevin noise added to both equations as in p5| ) are time-stepped 
using a standard explicit Euler method with a finite-difference representation of the spatial 
derivatives on a cubic lattice. In a second domain Z) 2 , consisting of the rest of the liquid, 
Eq. |5| reduces to the heat diffusion equation that is solved using a large ensemble of random 
walkers. The walkers make progressively larger steps with increasing distance away from the 
interface, and this adaptive stepping reduces the computer time in a way that is similar to 
adaptive meshing algorithms. The two methods are interfaced at the conversion boundary 
between domains D\ and D2 in such a way that heat is exactly conserved. 

All the details of this algorithm are described in Ref. |26j and need not be reproduced 
here. In addition, the finite-difference implementation of the phase- field equations in Di, 
which was used in Ref. [261, has been described in detail in Ref. [0. Here we make two 
additional improvements to this implementation. First, we integrate the equation for u (Eq. 
||D in D\ on a twice coarser mesh than the equation for (Eq. values of u on the fine 
mesh are then simply calculated by interpolating between values of u on the coarse mesh. 
Numerical tests show that this approach gives essentially identical results as using a single 
mesh. Second, we write the evolution equation for <fi derived from Eq. |] by inserting all the 
anisotropic factors in the form 



15 



? - w 2 



r (n) d t (f) = W$V 2 (j) + V ■ ([W(n 

+ (0-MA M )(l-0 2 ))(l-0 2 ) 



V0 



w=x,y,z 



(29) 



where the second and third terms on the right-hand side above vanish in the isotropic limit 
e c — > 0. The spatial derivatives associated with these terms are discretized using centered 



finite-difference expressions, as before ||22|| . In contrast, the isotropic Laplacian (first term 
on the right-hand-side above) is discretized using the 18-point formula 



1 



E^ + ^E^- 12 ^ 



(30) 



*j. fc 3Ax 2 

where J2nn 4> an d Unnn 4> denote the sums of 0- values evaluated at the 6 nearest-neighbor 
and the 12 next-nearest-neighbor positions on the cubic lattice, respectively, where Ax is 
the lattice spacing (equal along x, y, and z). The key property of this 18-point formula is 
that it makes the lattice corrections to the lowest order forms of the capillary and kinetic 
anisotropies defined by Eqs. [23| and ^ vanish at order Ax 2 . This can be seen by performing 
a Taylor expansion in Ax of the discretized Laplacian. For planar interfaces parallel to 
the (100), (110), and (111) planes this expansion reduces to the form V 2 ~ d 2 <f)/dn 2 + 
A(Ax) 2 <9 4 0/<9n 4 + 0(Ax 4 ), where n is the coordinate normal to the interface, and A is a 
numerical constant. For the 18-point formula, A is identical for all three planes. In contrast, 
A has a different value in each plane if the Laplacian is discretized using the standard finite- 
difference expression that only involves the nearest-neighbor sum. This, in turn, generates 
order Ax 2 lattice corrections to the capillary and kinetic anisotropies that can be taken 



into account following the estimation procedure developed in Ref. ||22| . The use of Eq. ^ 
eliminates the need for this procedure to resolve small physical anisotropies and pushes 
anisotropic effects due to the lattice discretization to a higher order (Ax 4 ) where numerical 
tests show that they are small enough to be neglected. 

Let us briefly comment on the sources of noise in our present numerical algorithm. 
The noise produced by the random walkers is only generated in D 2 , whereas the Gaussian 
thermal noise that is added through stochastic variables in the phase-field equations is only 
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generated in D\. It was shown in Ref. [£6J that the walker noise has a negligible effect on the 
interface evolution because its magnitude decreases by several orders of magnitude across the 
liquid layer surrounding the interface. Most importantly, this noise is insufficient to sustain 
sidebranching. One issue, however, is whether adding stochastic forces to the phase-field 
equations only in region Di (as opposed to the entire computation domain D 1 + D 2 , as would 
be the case if we had used a finite-difference scheme everywhere) reduces the strength of the 
thermal noise in a degree that affects the results. To investigate this issue numerically, we 
repeated certain test simulations with different values for the thickness of the liquid layer 
surrounding the interface, and found that the sidebranch amplitude remained approximately 
the same provided that this thickness is chosen large enough (about 25-50 times the width 
of the diffuse interface in the present study). Therefore, we conclude that our results are 
not significantly affected by only adding Langevin noise in D\. 

We use for material parameters of pure Ni: Tu = 1726 K, L = 2.311 x 10 9 J/m 3 , 
c p = 5.313 xlO 6 J/(m 3 K), D = 1(T 5 m 2 /sec, 7 = 0.326 J/m 2 , e c = 0.018, fi 100 = 0.52 m/sec, 
and /ino = 0.40 m/sec, where the capillary and kinetic parameters have been estimated by 
monitoring the fluctuations of the solid-liquid interface during MD simulations |3~iy4"7[] . From 
these values, we obtain L/c p = 435 K, do = 5.56 x 10~ 10 m, (3 = c p (l/^ 10 o + l//xno)/(2L) = 
5.084 x 10 -3 sec/m, e fc = (fi 100 - //no)/(^iio + A*iio) = 0.13, and F expt = 0.234. 

Simulations were carried out on a 800 x 520 x 520 lattice with the (100) growth direction 
along the larger dimension. We fully exploited the cubic symmetry of the crystal to reduce 
the computation time as described in f2~6fl . Furthermore, we used a spacing Ax = 0.8Wq 



and a time step At = 0.004 r , where W = 8AcV(5a/2) and r = 32A 2 d /V 25 - Note that 
the physical domain being simulated is actually the sum of the two aforementioned domains 
Di and D 2 , and D 2 extends farther than the lattice because the walkers make unrestricted 
off-lattice walks. We used as initial condition a small spherical seed and followed the mor- 
phological development of the interface up to the point where D\ is no longer contained 
inside the 3-d lattice. This allowed sufficient time for the dendrite tip to reach a steady- 
state velocity in the simulations for the parameters studied. Finally, we increased A from 2 
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to 16 as AT/(L/c p ) was decreased from 0.8 to 0.2, which yielded converged results that are 
reasonably independent of interface thickness (i.e. of A). 



IV. RESULTS AND DISCUSSION 
A. Comparison with experiments 

We compare in Fig. |2] the variation of the steady-state dendrite growth velocity with bulk 
undercooling obtained from the phase-field simulations with two sets of experimental mea- 
surements |]39|j42| . These measurements were made on small high-purity Ni droplets (a few 
millimeters in diameter) using a widely used experimental technique that consists of electro- 
magnetically levitating the droplets in order to avoid container-wall-induced heterogeneous 
nucleation, which is essential to attain very large undercoolings. 

In these experiments, crystallization of the undercooled melt is externally initiated by 
mechanical contact between the liquid droplet and a trigger needle at a well-defined under- 
cooling. After nucleation, the solidification front propagates outwards from the trigger point 
throughout the bulk of the droplet. The velocity of this front is measured by recording the 
abrupt change of surface temperature that is produced by its passage. Willnecker et al. 



measured the temperature with photodiodes in two regions of the droplet and obtained an 
average measure of the velocity from the time difference of the arrival of the front. Lum et 
al. m and Mat son [|4^] used a fast high-resolution video camera that images the thermal 
profile of the entire surface and that provides, in addition to the velocity, information about 
the morphological evolution of the solidification front. 

The comparison in Fig. || shows that the phase-field model predictions are in reason- 
ably good quantitative agreement with the measured values of the velocity up to a critical 
undercooling AT C m 200 K. For AT > AT C , the measured solidification velocity increases 
more slowly with undercooling than predicted by the model, especially for the data set of 



Lum et al. |41| and Matson [ 4^] . This break in the slope of the measured V — AT curve, 
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which is absent in the simulations, was first seen by Walker |38], and its existence has since 
been confirmed by numerous experiments [|39|-|42|[. 

From the analysis of the thermal imaging data, Lum et al. j|T] and Matson |4!| have 
concluded that this break is correlated with a change of morphology of the envelope of the 
solidification front, from angular for AT < AT C , to spherical for AT > AT C . Note that the 
spatial resolution of these images (64 x 64 pixels over an area of a few millimeter square) is 
far too coarse to image directly the interface on the scale of the dendrite tip, which is in the 
submicron range for these high undercoolings. Therefore, these images track the evolution 
of the macroscopic envelope of the solidification front, but not of the front itself. 

At these high growth rates, the thickness of the thermal diffusion boundary layer (which 
scales as D/V) is small enough for the tips of surviving secondary branches to grow inde- 
pendently, and at the same rate as primary branches. Moreover, since secondary branches 
grow at 90° from the trunks of primary dendrite branches sufficiently behind the tip, one 
would expect the envelope of the solidification front associated with steady-state dendrite 
growth along the [100] directions to be pyramidal, or angular once projected onto the droplet 



surface being imaged, as observed in Ref. ||42|| . In contrast, the observation of a transition 
to a smooth spherical envelope suggests that steady-state dendrite growth along a preferred 
crystallographic direction no longer exists for AT > AT C . 

Our phase-field simulations show no indication of the existence of such a transition with 
the present parameters used for Ni. The growth morphology remains dendritic, with a well- 
defined steady-state tip velocity, even for the largest simulated undercooling that is larger 
than AT C in the experiments. Examples of growth morphologies are shown in Figs. [|(a)-(b). 
Much larger size simulations, which are not feasible in three dimensions, would be necessary 
to describe the dynamics of the envelope of the solidification front on the scale where it has 
been imaged experimentally. However, we expect this envelope to be pyramidal over the 
entire range of undercoolings investigated here in view of the smooth V — AT curve and the 
persistence of a stable dendrite tip for AT > AT C . Larger size simulations conducted in two 
dimensions for the same parameters confirm this expectation. 
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B. Sensitivity to anisotropy and kinetic effects 



In order to investigate the source of the disagreement between simulations and exper- 
iments for AT > AT C , and also to gain insight into the factors controlling the interface 
dynamics, we carried out additional simulations by varying independently the magnitude of 
the capillary and kinetic anisotropy. Varying the magnitude of the capillary anisotropy (e c ) 
was found to have very little influence on the results. This is shown in Fig. 0] where we 



compare the phase-field simulation results for the value e c = 0.018 estimated for Ni [ TL\ and 
for e c = 0, corresponding to a purely isotropic interfacial energy. The dendrite velocities are 
almost identical over the whole range of undercooling investigated. 

In contrast, we found that both the solidification rate and the growth morphology depend 
sensitively on the magnitude of the kinetic anisotropy (e^) that was varied for a fixed value 
of the capillary anisotropy e c = 0.018. This is shown in Fig. |5] where we plot the computed 
dendrite velocity as a function of at a fixed undercooling. The velocity is seen to decrease 
as €k is lowered. Moreover, steady-state growth ceases to exist for less than a critical 
value el that varies from zero to about 0.02 for AT/(L/c p ) varying from about 0.2 to 0.7 as 
shown in Fig. [| We note that el should also depend on the capillary anisotropy that does 
not influence the interface dynamics when is large, as already mentioned, but that does 
compete with the kinetic anisotropy when is small (i.e. comparable to el). In particular, 
increasing e c shifts to higher undercoolings the region where steady-state dendrite growth 
ceases to exist in the e^-AT plane of Fig. ||. 

For 6^ < el, interface growth is dominated by tip-splitting, which produces new branches, 
and branch-termination by overgrowth from neighboring branches. Since the branches do 
not have any obvious preferred growth direction, the envelope of the interface on a scale 
larger than the average branch spacing becomes progressively more spherical as one moves 
away from the boundary that marks the lower limit of existence of dendrite growth in the 
efc-AT plane of Fig. |6|. This is illustrated in Fig. ^ that shows space-time plots of the 
interface for two extreme values of e& at fixed AT. An example of a 3-d spherical envelope 
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morphology is shown in Fig. 0(c). This structure is qualitatively similar to the so-called 
dense-branching morphology characteristic of diffusion-limited growth with tip-splitting JIB . 
This morphology has been observed in a wide range of systems, including two-dimensional 



viscous flow in a Hele-Shaw cell and amorphous annealing ||48|| . Fig. |] shows that the average 
velocity of the envelope of this morphology is much slower here than the tip velocity of a 
dendrite grown the same undercooling with a large kinetic anisotropy. 

The above results show that the interface dynamics is strongly influenced by kinetic 
effects. This is due to the fact that, for the large growth rates investigated here, the interface 
kinetic undercooling (AT/ = /iioo^O * s a significant fraction of the total bulk undercooling, 
whereas the curvature undercooling associated with capillary forces can be estimated to be 
one to two orders of magnitude smaller than ATj from the computed values of the tip radius. 
Of course, one would expect the opposite to be true for low undercoolings, but we have not 
investigated this range here for the parameters of Ni. 

C. Comparison with solvability theory and previous numerical studies 

The phase-field results are in good quantitative agreement with the predictions of an 
analytical solvability theory developed by Lee et al. that includes both capillary and kinetic 
effects [f49[,|50f . This theory is a direct extension of the linearized solvability theory developed 



previously by Barbieri and Langer [|14j], which included only capillary effects. In the theory of 
Lee et al, the steady-state growth equations are linearized around a paraboloid of revolution 
and solved using a WKB approximation. The resulting solvability condition for the existence 
of a smooth dendrite solution is expressed in the form of a vanishing integral condition where 
the integral is evaluated numerically. Details of this theory will be presented elsewhere and 
we only summarize here the main results. 

We compare in Fig. |] the dendrite velocity-undercooling curve predicted by this theory, 
for the same parameters of Ni as used in the phase-field model, to the simulation results. 
The agreement is remarkably good given the fact that the tip shape significantly deviates 
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from a paraboloid of revolution in this growth regime dominated by anisotropic kinetics. 
Furthermore, this theory predicts that dendrite growth ceases to exist in a region of the 
€k — AT plane that is below the solid curve in Fig. |6| and that is in reasonably good 
quantitative agreement with the phase-field results. The breakdown of the dendrite solution 
occurs because the main steady-state dendrite branch merges with a lower- velocity unstable 
branch at a critical undercooling. Hence, dendrites only exist below this critical undercooling 
as shown in the morphology diagram of Fig. || 

It should be emphasized that the transition from dendritic growth to dense-branching 
morphology with increasing undercooling is fundamentally different in the present kinetic- 
dominated regime than in the opposite capillary-dominated regime that has been previously 
investigated in both 2-d [[51]] (with earlier references therein) and 3-d |52]. In the latter 



regime, the disappearance of dendrites has been shown to be associated with noise-induced 
tip splitting [|51j, or with the appearance of other branches of steady-state solutions with a 
split tip |51|j5^]. In 2-d, the building block of the dense-branching or seaweed morphology 



is the so-called doublon consisting of two asymmetric fingers that grow cooperatively. 
In a capillary dominated growth regime, doublons exist above a critical undercooling that 
depends on the magnitude of the capillary anisotropy. Beyond that threshold, the main 
steady-state branch of dendrite solutions still exists, but doublons grow faster and outrun 
dendrites. In 3-d, for growth in a channel, a triplet structure consisting of three asymmetric 
fingers has been found in numerical simulations []52[ . No detailed quantitative study has 



been carried out, but it might be expected that, qualitatively, triplets play the same role in 
3-d as doublons in 2-d. In contrast, in the present kinetic-dominated regime, we have seen 
no evidence of doublons (triplets) in our 2-d (3-d) phase-field simulations. Consequently, 
the emergence of the dense-branching morphology is directly linked to the termination of 
the main dendrite steady-state branch at some critical undercooling |50[ . 



To further check this conclusion, we repeated without noise the phase-field simulations 
corresponding to isotropic kinetics and high undercoolings in the morphology diagram of Fig. 
H We found that tip splitting still occurred. Hence, we conclude that noise is not required 
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here to induce splitting. It is probable that noise shifts the critical undercooling at which 
splitting occurs, but we have not conducted an exhaustive (and numerically intensive) study 
to determine the magnitude of this shift. The relatively good agreement between solvability 
and phase-field results in Fig. ^| suggests that this shift is small. 

In a recent 3-d phase-field simulation study that focused on a capillary-dominated regime 
without kinetics ||55|| , we found that the dendrite tip was destroyed by thermal noise above a 
critical undercooling of about half the hypercooling limit that roughly coincides with AT C in 
the pure Ni experiments. Furthermore, we conjectured that this noise-induced tip splitting 
behavior could potentially explain the universally observed break in the V — AT curve, if 
it persists in the presence of kinetic effects 

The present study clearly demonstrates that noise-induced tip splitting is prevented by 
strongly anisotropic kinetics for the values of (3q and estimated in the MD simulations, 
which invalidates this hypothesis. It is possible that the values of these parameters could 
vary somewhat if the MD simulations were to be repeated with more accurate interatomic 
potentials. However, would have to be reduced by about one order of magnitude for a 
transition from dendrites to a spherical growth morphology to occur with increasing under- 
cooling according to the results of Fig. [| Moreover, since the phase-field results do agree 
quantitatively with experiments for AT < AT C , would have to decrease abruptly close to 
AT C to explain the data, which seems unlikely. 

Another possibility is that (3q is not correctly estimated by the atomistic simulations, 
and that the model could match the data over the whole undercooling range for a smaller 
value of (3q and a very weak kinetic anisotropy. This also seems unlikely because (3q has 
been extracted from MD simulations with different methods and potentials, which have all 



yielded results consistent with the Broughton-Gilmer- Jackson (BGK) model |55 |. What 
remains more uncertain is the magnitude of the kinetic anisotropy. 



Very recently, Mullis and Cochrane 53 carried out 2-d phase- field simulations of the 



crystallization of a pure melt using an estimate of (3q based on the Turnbull model [SB| that 



is several times smaller than the value predicted by MD simulations or the BGK model |35 
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They also used a very weak kinetic anisotropy induced by the orientation dependence of the 
interface thickness in the phase-field model. They found that, for these kinetic parameters, 
dendrites are destroyed by noise and replaced by a seaweed structure at large growth rate. 
Furthermore, they conjectured, as in Ref. ||55|| , that this noise- induced tip splitting could be 



the mechanism of the envelope transition observed by Matson [32 ] . Again, the present study 
invalidates this hypothesis, at least for pure melts, under the assumption that the kinetic 
parameters predicted by the MD simulations are correct. 

Finally, the disappearance of the dendrite branch was not seen in a previous numerical 



solvability calculation by Ben Amar []54j with isotropic kinetics (e& = 0). The reason is that 



a much smaller kinetic coefficient was used {(5qD /cLq ~ 15 in Ref. |54j] instead of (3qD/(1q ~ 90 



here), which extends the steady-state dendrite branch to higher undercoolings. 

D. Effects not included in the model 

If we assume that the existing MD estimates of interface kinetic parameters are correct, 
then the most interesting conclusion from the present study is that other physical effects 
neglected here influence the interface dynamics for AT > AT C . Since we have used a phase- 
field model that describes the crystallization of a pure substance with only diffusive transport 
in the melt, two possible candidates are fluid flow and the presence of very dilute impurities 
in the melt, which are both unavoidably present in experiments. 

An analytical study of flow has been carried out by Horvay [|56||, motivated by experi- 



mental observations by Walker J57J in pure Ni. He concluded that the negative pressure gen- 
erated by the density difference between solid and liquid should suffice to induce cavitation 
during the very initial stage of crystallization after nucleation. Although this phenomenon 
may be unrelated to the observed morphology transition, it raises the issue of whether fluid 
flow could modify sufficiently the conditions at the interface (during the subsequent stage 
where the entire melt crystallizes) to alter the MD estimates of kinetic parameters, and 
consequently the phase-field predictions. 
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Although the experiments to which we have compared our phase-field results used high 
purity Ni (e.g. better than 99.99% in Ref. |39|]), some small amount of contamination may 
play a role. Recently, Cochrane et al. have investigated the effect of the oxygen content on 



the velocity-undercooling curve in Cu [ 58| . They found that the typically observed break 
in the slope of this curve was absent for an oxygen content less than 200 ppm, but present 
for Cu doped with 600 ppm of O. Furthermore, there is significant variation in the V — AT 
curves reported from various experiments in Ni (see Fig. 1 in Lum et al), which could 
potentially result from oxygen contamination. 

Even very dilute amounts of impurities can strongly influence the interface dynamics 
because solutal diffusion in the liquid is about four orders of magnitude slower than heat 
diffusion. In addition, impurities are known to be trapped by the rapidly advancing interface 
for the velocities corresponding to the undercooling range studied here |59|. This solute 
trapping effect need not be isotropic. If it is sufficiently anisotropic, it could potentially alter 
the interface dynamics. To our knowledge, however, nothing is known either experimentally 
or theoretically about this anisotropy in metallic alloys. 

V. CONCLUSIONS 

The crystallization of highly undercooled Ni melts was modeled using a computationally 
efficient phase-field approach together with capillary and kinetic properties of the solid-liquid 



interface that have been recently estimated by atomistic simulations |3T1 , |32| . In particular, 
the value used here for the kinetic coefficient (flo) is much larger than the estimate of 
the Turbull model |36| that has been traditionally used to model dendrite growth at high 



undercooling. Consequently, the interface dynamics is much more strongly dominated by 
kinetics, and our results are drastically different from those obtained in previous studies 
based on the Turbull estimate [[41],[4^,^3|,|4j] or in related studies that focused on a purely 
capillary-dominated regime without kinetics pT| , |55*| . 

The dependence of dendrite growth velocity on undercooling predicted by the phase- 
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field simulations is in good quantitative agreement with experimental data for AT < AT C , 
where AT C is the critical undercooling corresponding to a break in the slope of the measured 
V — AT curve. For AT > AT C , the velocity increases more rapidly with undercooling in 
the simulations, which do not show this break, than in the experiments. Furthermore, the 
simulations predict that dendrite growth persists even for the largest undercoolings that 
have been studied experimentally. In contrast, the thermal imaging data of Lum et al. |41 



and Matson j|2| show that the envelope of the solid-liquid interface becomes spherical for 



AT > AT C , as opposed to angular for AT < AT C , which suggests that dendrites cease to 
exist for AT > AT C in apparent disagreement with the present results. 

Phase-field simulations in which the anisotropy parameters are systematically varied 
show that the growth rate and morphology depend sensitively on the kinetic anisotropy. In 
particular, they show that dendrites cease to exist above a critical undercooling comparable 
to AT C if the kinetic anisotropy is about one order of magnitude smaller than presently 
estimated by the atomistic simulations. For AT > AT C , these simulations yield a dense- 
branching morphology characteristic of tip-splitting-dominated growth. Furthermore, they 
show that the doublon (triplet) in 2-d (3-d) does not appear to be the underlying building 
block of this morphology in this strongly kinetic-dominated regime, in contrast to the sea- 
weed structure that has been simulated in 2-d without kinetics J5l| or with the much smaller 
value of (3q based on the Turnbull model (and noise) |53| . 



The phase-field results are in good quantitative agreement with the predictions of a 
linearized solvability theory that includes both capillary and kinetic effects ||50|| . This the- 
ory predicts dendrite velocities that are in remarkably good quantitative agreement with 
phase-field results for the large kinetic anisotropy value predicted by the MD simulations. 
Furthermore, it predicts that for a weak kinetic anisotropy, the steady-state dendrite branch 
terminates at some critical undercooling by merging with a lower velocity unstable branch. 
This provides a natural explanation for the disappearance of dendrites in the phase-field 
simulations with a sufficiently weak kinetic anisotropy, even without noise. 

These results emphasize the need for a precise knowledge of both the magnitude of the 
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kinetic coefficient {(3q) and its anisotropy (e^) to model reliably morphological evolution in 
a kinetic-dominated growth regime. There is still a significant uncertainty in the value of 
obtained by MD simulations |3l| . The present simulations, however, do show that must 
be large, as predicted JJl]], to reasonably fit the dendrite velocity data for AT < AT C . This 
is, of course, if one uses the value of /3 predicted by both MD simulations and the BGK 



model |35| , which is much less uncertain than the anisotropy. 

We are led to conclude that, if the present anisotropy value is confirmed, the widely 
observed velocity break in the V — AT curve [|38|-fH] is not reproduced by the standard 
model of the solidification of a pure melt with purely diffusive transport. This suggests that 
this break, and the associated envelope transition |[42|| , may be caused by flow effects or 
dilute impurities. The latter seems possible since even very dilute impurities in Cu have 



been shown to influence the dendrite velocity and grain structure |58|. An extension of the 
present study to link atomistic and phase-field simulations in alloys is currently in progress 
to test this hypothesis and its consequences for grain refinement [|60| . 
Acknowledgments 

This research was supported by the US Department of Energy through grant No DE- 
FG02-92ER45471 and additional funds provided by the Computational Materials Science 
Network. Computations were carried out at NERSC and NU-ASCC. We wish to thank Jeff 
Hoyt and Mark Asta for valuable discussions and input. 



27 



FIGURES 

FIG. 1. Plots of the dimensionless interface velocity v of an isothermal planar interface versus 
driving force — h in the phase- field model. 

FIG. 2. Dendrite velocity versus undercooling from two sets of experiments in Ni by Lum et al. 
|4l|| (open squares) and Willnecker et al. |3^] (open circles), and from the phase- field simulations 
(filled triangles and solid line) with A =16, 12, 8, 5, 4, 3, and 2 for AT/(L/c p ) = 0.2, 0.3, 0.4, 
0.5, 0.6, 0.7, and 0.8, respectively, and with the anisotropy parameters e c = 0.018 and = 0.13 
estimated by atomistic simulations for pure Ni [3lj32"|. 



FIG. 3. Snapshots of 3-d interfaces for different anisotropy parameters: (a) A = 30, = 0.13, 
e c = 0.018, and AT/(L/c p ) = 0.2; (b) A = 20, e k = 0.13, e c = 0.018, and AT/(L/c p ) = 0.7; (c) 
A = 20, efc = 0, e c = 0.018, and AT/(L/c p ) = 0.6. Larger values of A were used in these simulations 
to obtain more developed structures for illustrative purposes. 

FIG. 4. Dendrite velocity versus undercooling obtained from the phase-field simulations for 
e/c = 0.13 and e c = 0.018 (filled triangles and solid line), ej. = 0.13 and an isotropic interfacial 
energy e c = (open circles), and average velocity of the dense-branching morphology with 6^ = 
and e c = 0.018 (filled circle). The predictions of solvability theory (thick dashed line) for = 0.13 
and e c = 0.018 agree quantitatively with the simulations. 

FIG. 5. Dendrite tip velocity (filled triangles) and average velocity of the envelope of the 
interface (filled circles) versus kinetic anisotropy for A = 4, e c = 0.018, and AT/(L/c p ) = 0.7. 

FIG. 6. Growth morphology diagram showing the region of existence of dendrite (filled tri- 
angles) and dense-branching (filled circles) morphologies in the plane of kinetic anisotropy and 
undercooling for the phase-field simulations with e c = 0.018. For the same parameters, solvabil- 
ity theory predicts that steady-state dendrite growth solutions only exist above the solid line, in 
reasonably good quantitative agreement with the phase- field results. 
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FIG. 7. Plots of interface evolution for strongly anisotropic (a) and isotropic (b) kinetics. 
Interfaces in (a) and (b) were extracted from the same simulations that produced the dendrite and 
dense-branching morphology shown in Figs. ||(b) and|||(c), respectively. The interfaces represent 
cuts of the solid-liquid interface in a plane perpendicular to the [010] direction and are shown at 
equal time intervals. The frame width and height are both 8 fim in (a) and (b). 
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FIG.l: J.Bragard et al. 
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FIG. 2: J.Bragard et al. 
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FIG. 3a: J.Bragard et al. 
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FIG. 3b: J.Bragard et al. 
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FIG. 3c: J.Bragard et al. 
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FIG.4: J.Bragard et al. 




FIG. 5: J.Bragard et al. 
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FIG. 6: J.Bragard et al. 
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FIG. 7a: J.Bragard et al. 
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FIG. 7b: J.Bragard et al. 




43 



